lab1
1a Species Distribution Modeling- Exploratory Data Analysis
Install packages
# load packages, installing if missing
if (!require(librarian)){
install.packages("librarian")
library(librarian)
}
librarian::shelf(
dismo, dplyr, DT, ggplot2, here, htmltools, leaflet, mapview, purrr, raster, readr, rgbif, rgdal, rJava, sdmpredictors, sf, spocc, tidyr, geojsonio)
select <- dplyr::select # overwrite raster::select
# set random seed for reproducibility
set.seed(42)
# directory to store data
dir_data <- here("data/sdm")
dir.create(dir_data, showWarnings = F)Get species observations
American Black Bear
Ursus americanus (American Black Bear)
obs_csv <- file.path(dir_data, "obs.csv")
obs_geo <- file.path(dir_data, "obs.geojson")
# get species occurrence data from GBIF with coordinates
(res <- spocc::occ(
query = 'Ursus americanus',
from = 'gbif', has_coords = T,
limit = 10000))## Searched: gbif
## Occurrences - Found: 19,547, Returned: 10,000
## Search type: Scientific
## gbif: Ursus americanus (10000)(res_large_limit <- spocc::occ(
query = 'Ursus americanus',
from = 'gbif', has_coords = T,
limit = 100000))## Searched: gbif
## Occurrences - Found: 19,547, Returned: 19,547
## Search type: Scientific
## gbif: Ursus americanus (19547)# extract data frame from result
df <- res$gbif$data[[1]]
readr::write_csv(df, obs_csv)
anyDuplicated(df)## [1] 0# convert to points of observation from lon/lat columns in data frame
obs <- df %>%
filter(basisOfRecord == "HUMAN_OBSERVATION") %>%
sf::st_as_sf(
coords = c("longitude", "latitude"),
crs = st_crs(4326)) %>%
select(prov, key) # save space (joinable from obs_csv)
sf::write_sf(obs, obs_geo, delete_dsn=T)
obs <- sf::read_sf(obs_geo)
nrow(obs) # number of rows## [1] 9993# show points on map
mapview::mapview(obs, map.types = "Esri.WorldImagery")Question: How many observations total are in GBIF for your species?
There are 19547 total observations in GBIF for my species.
Question: Did you have to perform any data cleaning steps?
I do not see any odd observations in the data. The cleaning I had to do was minimal, and consisted only of filtering the basis of record to “HUMAN_OBSERVATION”
Get environmental data
dir_env <- file.path(dir_data, "env")
# set a default data directory
options(sdmpredictors_datadir = dir_env)
# choosing terrestrial
env_datasets <- sdmpredictors::list_datasets(terrestrial = TRUE, marine = FALSE)
# show table of datasets
#env_datasets %>%
# select(dataset_code, description, citation) %>%
#DT::datatable()# choose datasets for a vector
env_datasets_vec <- c("WorldClim", "ENVIREM")
# get layers
env_layers <- sdmpredictors::list_layers(env_datasets_vec)
#DT::datatable(env_layers)Question: What environmental layers did you choose as predictors? Can you find any support for these in the literature?
The layers I chose as predictors are WC_alt (altitude), WC_bio1 (annual mean temperature), WC_bio11 (mean temperature of coldest quarter), WC_bio10 (mean temperature of warmest quarter), WC_bio4 (temperature seasonality), ER_tri (terrain roughness index), ER_topoWet (topographic wetness), WC_bio12 (annual precipitation), WC_bio7 (annual temperature range), latitude, and longitude. Since my prior knowledge of black bear habitat is limited and comes mostly from common knowledge, I decided to choose a wide range of environmental factors as possible predictors.
# choose layers after some inspection and perhaps consulting literature
env_layers_vec <- c("WC_alt", "WC_bio1", "WC_bio11", "WC_bio10", "WC_bio4", "ER_tri", "ER_topoWet", "WC_bio12", "WC_bio7")
# get layers
env_stack <- load_layers(env_layers_vec)
# interactive plot layers, hiding all but first (select others)
# mapview(env_stack, hide = T) # makes the html too big for Github
plot(env_stack, nc=2)
obs_hull_geo <- file.path(dir_data, "obs_hull.geojson")
env_stack_grd <- file.path(dir_data, "env_stack.grd")
# make convex hull around points of observation
obs_hull <- sf::st_convex_hull(st_union(obs))
# save obs hull
write_sf(obs_hull, obs_hull_geo)
obs_hull <- read_sf(obs_hull_geo)
# show points on map
mapview(
list(obs, obs_hull))obs_hull_sp <- sf::as_Spatial(obs_hull)
env_stack <- raster::mask(env_stack, obs_hull_sp) %>%
raster::crop(extent(obs_hull_sp))
writeRaster(env_stack, env_stack_grd, overwrite=T)
env_stack <- stack(env_stack_grd)
# show map
# mapview(obs) +
# mapview(env_stack, hide = T) # makes html too big for Github
plot(env_stack, nc=2)
## Pseudo-Absence
absence_geo <- file.path(dir_data, "absence.geojson")
pts_geo <- file.path(dir_data, "pts.geojson")
pts_env_csv <- file.path(dir_data, "pts_env.csv")
# get raster count of observations
r_obs <- rasterize(
sf::as_Spatial(obs), env_stack[[1]], field=1, fun='count')
# show map
#mapview(obs) +
mapview(r_obs)# create mask for
r_mask <- mask(env_stack[[1]] > -Inf, r_obs, inverse=T)
# generate random points inside mask
absence <- dismo::randomPoints(r_mask, nrow(obs)) %>%
as_tibble() %>%
st_as_sf(coords = c("x", "y"), crs = 4326)
write_sf(absence, absence_geo, delete_dsn=T)
absence <- read_sf(absence_geo)
# show map of presence, ie obs, and absence
mapview(obs, col.regions = "green") +
mapview(absence, col.regions = "gray")# combine presence and absence into single set of labeled points
pts <- rbind(
obs %>%
mutate(
present = 1) %>%
select(present, key),
absence %>%
mutate(
present = 0,
key = NA)) %>%
mutate(
ID = 1:n()) %>%
relocate(ID)
write_sf(pts, pts_geo, delete_dsn=T)
# extract raster values for points
pts_env <- raster::extract(env_stack, as_Spatial(pts), df=TRUE) %>%
tibble() %>%
# join present and geometry columns to raster value results for points
left_join(
pts %>%
select(ID, present),
by = "ID") %>%
relocate(present, .after = ID) %>%
# extract lon, lat as single columns
mutate(
#present = factor(present),
lon = st_coordinates(geometry)[,1],
lat = st_coordinates(geometry)[,2]) %>%
select(-geometry)
write_csv(pts_env, pts_env_csv)
pts_env <- read_csv(pts_env_csv)
pts_env %>%
# show first 10 presence, last 10 absence
slice(c(1:10, (nrow(pts_env)-9):nrow(pts_env))) %>%
DT::datatable(
rownames = F,
options = list(
dom = "t",
pageLength = 20))Term plots
pts_env %>%
select(-ID) %>%
mutate(
present = factor(present)) %>%
pivot_longer(-present) %>%
ggplot() +
geom_density(aes(x = value, fill = present)) +
scale_fill_manual(values = alpha(c("gray", "green"), 0.5)) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
theme_bw() +
facet_wrap(~name, scales = "free") +
theme(
legend.position = c(1, 0),
legend.justification = c(1, 0))
1b. Logistic Regression
librarian::shelf(
DT, dplyr, dismo, GGally, here, readr, tidyr)
select <- dplyr::select # overwrite raster::select
options(readr.show_col_types = F)
dir_data <- here("data/sdm")
pts_env_csv <- file.path(dir_data, "pts_env.csv")
pts_env <- read_csv(pts_env_csv)
nrow(pts_env)## [1] 19986datatable(pts_env, rownames = F)